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Abstract: The virtual two-loop corrections for Higgs production in gluon fusion are 
calculated analytically in QCD for arbitrary Higgs and quark masses. Both scalar and 
pseudo-scalar Higgs bosons are considered. The results are obtained by expanding the 
known one-dimensional integral representation in terms of mn/mq, and matching it with 
a suitably chosen ansatz of Harmonic Poly logarithms. This ansatz is motivated by the 
known analytic result for the Higgs decay rate into two photons. The method also allows 
us to check this result and to extend it to the pseudo-scalar decay rate. 
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1. Introduction 

The gluon fusion process for Higgs production at a hadron collider has been studied in 
great detail over the last few years (for a recent review, see Ref. [1]). It is well-known to 
be the dominant mode in the Standard Model and also in most of the usually considered 
supersymmetric parameter space. The fact that the next-to-leading order QCD correc- 
tions [2-4] increase the cross section by more than 70% triggered more detailed studies 
of higher order effects. In particular, the NNLO [5-9] and quite recently even the leading 
threshold-enhanced N 3 LO [10] corrections were evaluated in the heavy-top limit, indicating 
a well-behaved perturbative expansion of the total cross section. Meanwhile, the NNLO 
effects are known also for differential quantities in terms of a partonic NNLO Monte Carlo 
program [11,12], allowing to simulate experimental cuts, for example. 

In contrast to the NNLO calculations which currently all rely on the heavy-top limit, the 
inclusive NLO effects were calculated for arbitrary values of the Higgs boson mass and the 
mass of the quark that mediates the gluon-Higgs coupling [4, 13, 14]. In fact, it is this 
calculation that justifies the use of the heavy-top limit at NNLO, because it explicitely 
demonstrates the excellent quality of this limit even at Higgs masses close to the quark 
threshold run ~ 2m q and beyond. Probably the most important application of the general 
mn/mq dependence currently is super symmetry, where bottom quarks can contribute sig- 
nificantly to the gluon-Higgs coupling due to a potential enhancement proportional to tan (3 
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of their Yukawa coupling to Higgs bosons. For bottom quarks, an analogous "heavy-quark" 
approximation would certainly be very doubtful in this context [4, 14-16]. 

Considering the importance of the full mass dependence, it is somewhat surprising that the 
status of the NLO calculation is still at the level of more than ten years ago. By then, the 
result was obtained in terms of a rather lengthy one-dimensional integral representation 
and implemented in a FORTRAN routine. On the one hand, this makes it rather difficult 
to import the result into other programs, of course. On the other hand, it is practically 
impossible to further manipulate the result. 

The lack of an analytical result is also surprising in view of the great technical progress 
since the original work of Ref. [4]. In fact, the corresponding Feynman integrals belong 
to a class that currently receives great attention due to its importance for electro-weak 
precision observables (see, e.g., Refs. [17-19] and references therein). It turns out indeed 
that all integrals needed for a representation of the 2-loop virtual terms in closed form have 
been evaluated in the literature. 

In this paper we derive this analytic formula. Let us stress though that we did not evaluate 
the corresponding Feynman diagrams; rather, we used the integral representation given in 
Ref. [4] and evaluated it analytically. The method we followed is rather unconventional 
but not new. For the sake of brevity, we will refer to it as Expansion and Inversion (E&l) 
in what follows. It relies on the identity theorem for power series: Two analytic functions 
are the same if their Taylor series are the same. A more detailed description of the method 
and its realization will be given in Section |2[ 

2. Discussion of the method; calculation of the decay rates 

The idea behind our approach is that, if two physical processes correspond to a similar set of 
Feynman diagrams (kinematics, mass assignment), their cross sections should be described 
by a common set of analytical functions. Thus, if one processes is known, one can establish 
an ansatz for the other one by a linear combination of these functions, with unknown 
coefficients. In the E&I method, one then evaluates the power series of the unknown cross 
section in a certain limit and compares it with the corresponding expansion of the ansatz. 
This leads to a system of linear equations for the unknown coefficients which can be solved 
uniquely if the depth of the expansion matches the number of unknowns. In general, it is 
advisable to overdetermine the system in order to confirm that the ansatz is complete. The 
identity theorem for power series ensures that the solution obtained in this way is indeed 
the analytical result for the cross section. 

It is important to realize that, while the intermediate power series approximates the full 
result only within the radius of convergence, the final result is valid for arbitrary values of 
the parameters. Thus, the comparison of the final result with a numerical evaluation of 
the original integral outside the radius of convergence provides one of the most powerful 
checks on the calculation. 
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The main advantage of the E&I approach is that, in most cases, the power series of the 
cross section can be obtained in a rather simple manner. A powerful tool for this goal is 
provided by asymptotic expansions of Feynman diagrams (see Refs. [20-22] and references 
therein). This method works directly at the level of Feynman integrals and has been fully 
automated for the case of Euclidean external momenta [23,24]. 

Very often, however, one can derive one-dimensional integral representations over finite 
integration regions. This can be achieved by introducing Feynman parameters, for example, 
and performing all but one integral analytically. If the interchange of integration and 
differentiation is possible, the power series expansion can be performed directly on the 
integrand, and the resulting integrals are in general much simpler than the original ones. 

Another example where the integrations are over finite regions is given by phase space 
integrals, and in fact, the E&I method has been used for the evaluation of the three- 
particle phase space integration in the case of Higgs production and the Drell-Yan process, 
both at NNLO [5,7,25,26]. 

In this paper, we apply the E&I method to obtain analytic formulae for the NLO predictions 
of the Higgs decay rate into two photons, as well as to the virtual two- loop corrections for 
Higgs production in gluon fusion. We consider both scalar and pseudo-scalar Higgs bosons 
such that our results are relevant also in supersymmetric scenarios or other extensions of 
the SM. 

For all these quantities, a one-dimensional integral representation is known [4]. By in- 
terchanging differentiation and integration, we can derive their power series in terms of 
mn/nT-q- A closed analytical result is only known for the NLO decay rate of a scalar Higgs 
boson into photons [27] . We use it as the main motivation of our ansatz in order to derive 
closed analytical expressions for the other quantities as well. It will be useful for the rest 
of this paper to quote the explicit result at this point. 



2.1 Decay rate H — ► 77 

The decay rate of a Higgs boson into two photons through NLO can be written as (see, 
e.g., Ref. [4]) 



GFa 2 m 3 H 



T(H - 77) 



128V2V 3 



Q 2 iA?{n) + 3 £ Q 2 q Af(r q ) + A»{tw) 



(2.1) 



where Gf is the Fermi constant, a is the electromagnetic fine-structure constant, mu 
denotes the Higgs mass, and Q q j the electric charge of quark q and lepton / in units of the 
proton charge. The variables Tj are defined as 



m 



H 



4m? ' 



(2.2) 



where rrii denotes the mass of particle i. Here and in what follows, m q = m q (n) denotes 
the MS quark mass renormalized at a mass scale ji. 
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H/A 



(a) (b) 
Figure 1: Feynman diagrams contributing to the decay rate for H/A — > 77 at LO. 



The amplitudes Af (r) and ^4^(r) arise from closed lepton and VF-boson loops, respectively 
(cf. Fig. [l]), and do not receive QCD corrections. They are given by [28,29] 



A?(r) = ^[r + (r-l)f(r)] =^Ff(r), 



1 



[2r 2 + 3r + 3(2r-2)/(r)], 



(2.3) 



where 
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(2.4) 




(a) (b) 
Figure 2: Sample Feynman diagrams contributing to the decay rate for H/A — > 77 at NLO. 

Ag (t) originates from the quark mediated photon-Higgs coupling, Fig. 0(a) and Fig. [2]. 
Through NLO QCD, one can write it as 

,2- 



A?(r) = iFf(r) 1 + ^ ( Cf(r) + Of (r)ln ^ 



7T 
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(2.5) 



with _Fq*(t) from Eq. ( |2.3| ). follows directly from the NLO renormalization group equa- 
tion 



It reads 



^ C? = ^ [t + (r - 2) /(r) - (r - l)r /'(r)] 



(2.6) 



(2.7) 



The analytical expression for has been obtained in Ref. [27] : x 

1 Eqs. (10) and (12) of Ref. [27] contain typos; thanks to O. Tarasov for immediate confirmation. Note 
that in order to compare Eq. J2.S{ ) with the formula in Ref. [27], one needs to use the identity \As{8 2 ) = 
4[Li 3 (0)+Li 3 (-0)]. 
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(2.10) 



s/i - T^ 1 + 1 ' 

For analytic continuation, it is always understood that r — > t + zO. 

The products of logarithms and polylogarithms in Eq. (|2.8|) can be expressed in terms 



of Harmonic Polylogarithms [30] of the form H(n;0), where n is an n-tuple with entries 
("indices") ±1 or 0. One finds that n < 4, and that at most one index is different from 
zero. 

This suggests to construct our ansatz from Harmonic Polylogarithms of this form, multi- 
plied by rational functions 



Pn{9) 



(1 



(2.11) 



where P n {9) is a polynomial in of degree n with unknown coefficients. In order to solve 
the resulting system of linear equations, one has to adjust the integer parameters n, k such 
that a suitable balance is obtained between the universality of the ansatz and the depth of 
the power series expansion that is required to determine the unknown coefficients. 



As a warm-up, we may try to reproduce Eq. ( |2.8[) from the one-dimensional integral rep- 
resentation of Ref. [4] using our approach. To this aim, we expand the integrands of 



- 5 - 



defined in Eqs. (A. 9) to (A. 13) of Ref. [4], 2 around the limit r = 0, keeping 
terms through order r 100 . 

It is clear that due to the complexity of the integrands and the required depth of the 
expansion we need to use efficient computer algebra tools. We found that the Taylor 
package [31] for Reduce [32] is particularly well suited for this kind of operations. In most 
cases, the results obtained from Taylor were checked against our own implementation of 
the relevant power series in Form [33]. The capabilities of Mathematica [34], on the 
other hand, are clearly not suited for expanding expressions of this complexity. 3 

For illustration of the method, let us consider the simplest one of the relevant integrals: 

' 1 — px(l — x) \ 



(2.12) 



where 
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Expanding the integrand around r = leads to very simple integrations in x, 
•l ( 



h 



dx { 2 x In x + t 



8x 2 - 8x 3 + lnx(l0x 2 - -x J 



56 x 3 - 



248 



304 x 4 - 



3 

1744 



x H x 5 + In x I 



136 



68 



32 



x° + 



5504 



15 



v 3 x ~Y X + T X " 

448 7 

x 



, /592 4 2128 

+ In x\ x 

V 3 15 



= 1184 fi 128 7 
x H — x — x 



15 



4480 



x 5 - 



156256 



45 



•x u + 



16192 



x 7 - 



7 

164704 
105 



■x 8 + 



97408 
315 



, / 12608 . 3904 fi 67712 7 
+ m x x x H — x — 



V 15 



105 



2080 x 8 + 512 ^ 
7 9 j 



+ 



(2.13) 



such that 
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The remaining integrals in Ref. [4] are more complex, but the integration of their expansion 
in r can always be evaluated in an elementary way. 



2 Thanks to M. Spira for clarification concerning some typos in the formulas of Ref. [4]. 
3 We remark that Mathematica 5.1 even produces a wrong result when expanding Li2(l — x) around 
= (also for Li3 etc.); a bug report has been submitted and acknowledged. 
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We note in passing that for the coefficient of r 100 , as it is required by our general ansatz, 
the integers in the numerator and the denominator are roughly of order 10 180 ; this should 
give an idea of the intermediate expressions' complexity. 

Nevertheless, this expression, together with the corresponding expansion of the ansatz, can 
be fed into Mathematica in order to solve the resulting system of linear equations. One 
finds 
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Needless to say that the result obtained for in this way is in agreement with Eq. fl2.g|), 
thus proving the consistency of the analytical result of Ref. [27] and the integral represen- 
tation of Ref. [4] . 

Let us add a few more remarks concerning the construction of the system of linear equa- 
tions. It happens that the structure of the ansatz can be restricted already from general 
considerations: The expansions of the entities calculated in this paper all consist only of 
integer powers of r multiplied by rational coefficients. The expansions of Harmonic Poly- 
logarithms of argument 9, on the other hand, contain noninteger powers of r, irrational 
numbers like £„, logarithms of r, and have a non- vanishing imaginary part. The fact that 
such terms do not appear in the expansions of the integrals to be matched produces a lot 
of equations that constrain our ansatz, regardless of the specific integral to be calculated. 



2.2 Decay rate A — > 77 



After confirming the analytical result for T(H — * 77) through NLO, we are now ready 
to apply the method to the pseudo-scalar case. Assuming the Minimal Supersymmetric 
Standard Model (MSSM) as underlying theory, we write, in analogy to Eq. (|2.1|): 
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32^7r 3 

with the lepton (I) and chargino (x^ 1 ) induced amplitudes 



, (2.16) 



Af(r) = Aj ± (r) = ^ ^ F A (r) 



T 



(2.17) 



where f(r) has been defined in Eq. ( |2.4| ). Note that there is no contribution from the W 
as loop particle due to CP invariance. In Eq. ( [yj) , vtla is the mass of the pseudo-scalar 
Higgs, and the r- variables are defined according to Eq. (|2.2|), with uia instead of ran. The 



specific values of the couplings gf- - ± are irrelevant for our analysis; they can be found in 
Ref. [4]. 
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In analogy to Eq. ( |2,5[) , the quark-induced amplitude is written as 



K{r) = Fq(t) 



IT 



m , 



(2.18) 



where again C 2 can be derived through a renormalization group equation analogous to 
Eq. 0: 



F A C 2 A = ^[f(r)-rf'(r)] . 



(2.19) 



Remarkably, when using the same ansatz as in the scalar case of Sect. 2.1, the resulting 



system of linear equations has no solution. The necessary generalization is to allow for 
terms ~ (1 + 9)~ 1 multiplying the HPLs, reflecting the well-known threshold singularity in 
the pseudo-scalar case at txia = 2m q . Once this is done, the E&I approach yields 
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3. Virtual corrections for grgr — >■ H/A 



An interesting application of our method is the analytical evaluation of the virtual two-loop 
corrections for Higgs production in gluon fusion for arbitrary values of the quark and Higgs 
boson mass. 
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Following Ref. [4], we write the inclusive NLO cross section as 

a(pp -» $ + X) = <j* 

<De{F,,4}, 

where r$ = m\/s with the center-of-mass energy s, and 
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with the gluon density functions g{x,Hp), depending on the factorization scale fip. The 
normalization factors are 
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(3.3) 



with Fq^ A defined in Eqs. (|2.3| ), ( 2.17| ). denotes the contributions from the virtual 
two-loop corrections, regularized by the infrared singular part of the cross section for real 
gluon emission (see Ref. [4] for details). It can be decomposed into 



C*=7r 2 + c* + 2/3 ln-^, 

mi 



(3.4) 



where (3$ = 11/4 — nj/6 is the lowest-order (3 function of QCD for rif active quark flavors. 
The Act* denote the contributions from radiation of quarks and gluons with initial state 
partons i,j G {q,q, <?}• At NLO perturbation theory, they correspond to massive one- 
loop three- and four-point functions which can be evaluated analytically using standard 
techniques [35] (see also Ref. [36]). They shall not be considered any further in this paper. 

The coefficient c* of the virtual corrections in Eq. ([^J) is parameterized as 
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(3.5) 



Similar to the decay rates, Bf follows from renormalization group considerations: 

B*(r) = 2C*(r), <5>£{H,A}, 



(3.6) 



with Cf from Eqs. (2/7) and ( [2,19|) . The factor of 2 arises from the fact that C* in Eq. ( |3.1| ) 
is defined at the level of the cross section rather than the amplitude. 

Bf is known again in terms of one-dimensional integrals, filling several pages (ii, . . . ,ig 
in App.A,B of Ref. [4]). Following the method described in Sect.||, we expanded the 
integrands for Bf(r) around r = up to order r 100 and mapped them onto a set of 
suitably chosen basis functions. 
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Figure 3: (a) Feynman integral contributing to the production rate gg 
decay rate H/A — > 77; it arises due to the self-coupling of gluons, see (b). 



H/A but not to the 



In the case of gluonic Higgs production, a new class of integrals occurs that cannot be 
expressed in terms of the functions used for the decay rates. The corresponding scalar 
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diagram is shown in Fig. ^ (a); it arises from the physical process due to the self- interaction 
of gluons, see Fig. |3](b). The scalar integral has been evaluated in Ref. [37]. The result 
contains a Harmonic Polylogarithm of weight four with two indices different from zero. We 
therefore enlarge our basis to include this kind of functions. Apart from that, the method 
works exactly like in the case of the decay rates, described in Sect.^. 

With defined in Eq. (glcD , we find for the scalar case 
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For the pseudo-scalar case we get 
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One may notice that the Harmonic Polylogarithm appearing in Ref. [37] is different from 
the one contributing to and B^. However, this is only due to an arbitrariness when 
choosing a basis of Harmonic Poly logarithms. In fact, 
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as can be seen using Appendix B of Ref. [38], for example. Using H(l, 0, —1, 0; 9), all S 2)2 
and Si 2 cancel in our final result. 



4. Numerical Results 



As already mentioned in Sect.||, the comparison of the final result with a numerical evalu- 
ation of the original integral provides one of the most important checks of the calculation. 

The solid and dashed lines of Fig. || (a) and (b) show the real and imaginary part of Fq 
and CiFq, respectively. The dotted lines in Fig. ^ show the results obtained from the 
intermediate power series when keeping terms of order r n with n = 10, 30, 90. Note that, as 
expected, the power series does not converge towards the analytic result beyond the radius 
of convergence, given by r = 1. In particular, the imaginary part is always zero. Thus, 
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Figure 4: Two-loop contributions to the (a) scalar and (b) pseudo-scalar Higgs decay rate into 
photons. The solid and dashed lines show the real and the imaginary part, respectively. The dotted 
lines correspond to the power series expansion up to order r™ with n = 10, 30, 90. 



for r > 1, the result arises solely from analytic continuation of the terms reconstructed 
through E&I. 

In addition, we were able to reproduce Figs. 5 and 18 of Ref. [4], using our results, to 
perfect agreement. We find a similar picture for the virtual corrections to gluon fusion, 
shown in Fig.||(a) and (b) for the scalar and the pseudo-scalar case, respectively. The 
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Figure 5: Infrared regularized virtual two- loop corrections to the (a) scalar and (b) pseudo-scalar 
Higgs production rate through gluon fusion. The solid and dashed lines show the real and the 
imaginary part, respectively. The dotted lines correspond to the power series expansion up to order 
t™ with n = 10, 30, 90. 



numerical evaluation of H(1,0, —1,0;$) was done with the help of the Mathematica file 
HPL4.m in [39]. 
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5. Conclusions 



The two-loop QCD results for the decay rate of a scalar or pseudo-scalar Higgs boson 
into photons, H/A — ► 77, as well as for the virtual corrections to the production modes 
gg — ► were presented in closed analytical form. In order to obtain these result, we 

first expanded the known one-dimensional integral representations in terms of small Higgs 
masses, and subsequently mapped this expansion onto a set of analytic functions. The final 
results, both for their real and imaginary part, are valid for arbitrary values of the quark 
and Higgs boson mass. They contain only polylogarithms or simpler functions and, in the 
case of gg — > H/A, one Harmonic Poly logarithm. 

Our formulas should be useful for implementations into physics analysis programs, or for 
quickly obtaining analytical limits to arbitrary accuracy. 

Let us finish by pointing out that the E&I method, in various flavors, has been quite useful 
already in the past (see Refs. [7,25,40-43], for example). Its combination with asymp- 
totic expansions may even carry the potential for an algorithmic evaluation of Feynman 
integrals. However, this not only requires much more efficient computer algebra tools for 
the expansion of Feynman diagrams. The more important task is to find suitable bases 
for certain classes of Feynman integrals. We believe that this is certainly a task worth 
pursuing. 
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